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Abstract. Using a large scale molecular dynamics computer simulation we investigate 
the dynamics of a supercooled melt of SiC-2 . We find that with increasing temperature 
the temperature dependence of the diffusion constants crosses over from an Arrhenius- 
law, with activation energies close to the experimental values, to a power-law depen- 
dence. We show that this crossover is related to the fact that at low temperatures the 
dynamics of the ions is dominated by hopping processes, whereas at high temperatures 
it shows the continuous flow-like motion proposed by the ideal version of mode-coupling 
theory (MCT). Finally we show that at low temperatures the dynamics of the system 
in the (3— relaxation regime obeys the factorization property, in agreement with MCT. 



INTRODUCTION 

In the last few years a significant advance in our understanding of the dynamics 
of glass-forming liquids has occurred. By now it is well established that features 
like, e.g., the cage-effect, the slow dynamics in the j3— relaxation regime, and the 
stretching of the relaxation curves in the a— relaxation regime, are characteristic for 
the dynamics of supercooled liquids, and that many of these features are described 
well, or have even been predicted, by mode-coupling theory (MCT) [1]. For fragile 
glass-formers [2] there are even systems for which MCT describes the relaxation 
dynamics correctly not only on a qualitative level, but even on a quantitative one, 
such as colloidal particles or Lennard- Jones systems [3-5]. The dynamics of strong 
glass- formers is, however, understood in much less detail since it is, e.g., not yet 
clear whether the features predicted by MCT to be present in the (3— relaxation 
regime, such as the critical decay or the von Schweidler law, do not exist in these 
system or are just obscured by other dynamical features, such as the boson-peak. 
In the present paper we therefore present some results of a large scale computer 
simulation in which we investigated the dynamics of a silica melt, the prototype of 
a strong glass former. 



MODEL AND DETAILS OF THE SIMULATION 



The silica model we study has been proposed by van Beest et al. (BKS) [6] and 
has been found in the past to reproduce well various structural properties of real 
silica [7]. In this potential the interaction between the ions is given by 

= ^ + Aj eM-Bijn,) - ^ . (i) 

r n 'ij 

Here is the distance between ions % and j, and the values of the partial charges 
qi and the constants Aij, Bij, and CV,- can be found in Refs. [6,7]. The simulations 
have been done at constant volume using; 8016 ions and a box size of 48.37A, thus 
at a density of 2.37g/cm 3 , close to the experimental value of the density. (Such a 
large system was needed to avoid finite size effect in the dynamics [8].) Using Ewald 
sums to evaluate the Coulombic part of the potential, the equations of motion have 
been integrated with the velocity form of the Verlet algorithm with a time step 
of 1.6 fs. The temperatures investigated were 6100 K, 5200 K, 4700 K, 4300 K, 
4000 K, 3760 K, 3580 K, 3400 K, 3250 K, 3100 K, 3000 K, 2900 K and 2750 K, and 
in order to improve the statistics of the results we averaged at each temperature 
over two independent runs. At each temperature the system was first equilibrated 
over a time span which significantly exceeded the typical relaxation time of the 
system at this temperature. At the lowest temperature the runs had 13 million 
time steps, which is equivalent to about 20 ns of real time. The total CPU time 
used for this temperature was 13 years of single processor time on a CRAY-T3E. 

RESULTS 

The simplest quantity which characterizes the dynamics of a system is (r 2 (t)), 
the mean squared displacement (MSD) of a tagged particle, which is calculated 
from the positions fj(t) by means of 

(r 2 (t)) = —Y,(\m-m\ 2 ) with a G {Si, O}. (2) 

a j=i 

In Fig. 1 we show the time dependence of (r 2 (t)} for the oxygen atoms for all 
temperatures investigated. (Qualitatively similar curves are found for the silicon 
atoms [9].) From this figure we see that at high temperatures the MSD crosses 
over from the ballistic motion of the particles, (r 2 (t)) oc t 2 , to the diffusive motion, 
(r 2 (t)) oc t. At low temperatures we see in addition to these two regimes also a third 
one at intermediate time scales, in which the MSD stays essentially constant over 
2-3 decades in time. In this time regime the particle is trapped in the cage formed 
by the atoms that surround it, and it takes the particle a long time before it is able 
to leave this cage. One of the main features of MCT is to give a self-consistent 
description of this motion, and below we will discuss some of these predictions. 
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FIGURE 1. Time dependence of the mean squared displacement for the oxygen atoms for all 
temperatures investigated. 

(Self-consistency is needed since the particles forming the cage are trapped by their 
neighbors.) 

Using the Einstein relation D = lim t _ >00 (r 2 (t))/6t, it is simple to calculate the 
diffusion constant D from the MSD. The temperature dependence of D is shown 
in Fig. 2 for and Dq in an Arrhenius plot. We see that at low temperatures 
the diffusion constants show the Arrhenius dependence expected for a strong glass- 
former. The activation energies found, 5.18 eV and 4.66 eV for silicon and oxygen, 
respectively, are in very good agreement with the experimental values of Brebec 
et al. (silicon: 6 eV) and Mikkelsen (oxygen: 4.7 eV) [10]. Since a similar good 
agreement has been found for the temperature dependence of the viscosity [11] we 
thus conclude that the present model of silica is not only able to reproduce reliably 
the static properties of amorphous silica [7], but also the dynamic ones. 

From the figure we also see that at high temperatures significant deviations 
from the Arrhenius law are observed. Such deviations are often found in fragile 
glass-formers and have been explained by means of MCT, which predicts that in 
supercooled systems there exists a critical temperature T c in the vicinity of which 
transport coefficients, such as the diffusion constant or the a— relaxation time r, 
show a power-law dependence on temperature, i.e. 
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FIGURE 2. Arrhenius plot of the diffusion constants. The solid lines are Arrhenius fits to our 
data at low temperatures. The dashed lines are power-law fits to our data at high temperatures 
with an exponent 7. 



where the exponent 7 is a system universal parameter, i.e. is independent of the 
transport coefficient. Motivated by this theoretical prediction we thus have at- 
tempted to fit the temperature dependence of the diffusion constants with such 
a power-law and have included the results of these fits in Fig. 2 as well (dashed 
lines). Note that in these fits the critical temperature T c was not used as a fit pa- 
rameter but fixed to 3330 K, a value we obtained from fits to the relaxation times 
t of the intermediate scattering function for different values of the wave-vector q 
(which will be discussed below). From the figure we can see that the so obtained 
fits reproduce the data very well over a temperature range in which the diffusion 
constants change by about 1.5 decades. This range is significantly smaller than 
the one typically found in fragile glass-formers (3-4 decades [12]) and it is often 
argued that such a small range is related to the dominant presence of activated 
processes, the so-called hopping processes, which are taken into account only in the 
"extended" version of MCT [13] and not in the "ideal" version [1]. 

One way to see the mentioned hopping processes directly is to study the self part 
of the van Hove correlation function, G s (r, t), which is given by 

1 N a 

G:(r,t) = —j:(5(r-\m-f t (0)\)} « G {Si, 0} . (4) 




FIGURE 3. r— dependence of the self part of the van Hove correlation function for the oxygen 
atoms. The times are 0.16 ps, 1.44 ps and then increase by a factor of about 2.3. The last curves 
are 15.2 ns and 19.2 ns . Inset, the same curves in a lin-log plot. The arrow marks the location 
of the first peak in the radial distribution function between two oxygen atoms. 

Thus Anr 2 G^ (r, t) is the probability to find a particle at time t at a distance 
r away from the place it was at t — 0. In Fig. 3 we show the r— dependence of 
A-nr 2 G° (r, t) for different times for the lowest temperature investigated. 

In the main figure we see that at intermediate times this function is peaked at 
around 0.4A and hardly depends on time. This observation can be understood 
by recalling that in this time regime, the f3— regime, the particles are trapped 
in the above discussed cages formed by the surrounding particles. Only for times 
corresponding to the a— relaxation this peak starts to decrease significantly and the 
function begins to show a secondary peak, or rather a shoulder, at a distance around 
2.6A, the location of the first peak in the radial distribution function between two 
oxygen atoms [7,14]. This secondary peak can be better seen in a semi-logarithmic 
plot of this function, which is shown as an inset of the figure. Hence we conclude 
that at low temperatures the cage for the oxygen atoms breaks up because an atom 
jumps to the location of one of its six nearest oxygen neighbors, assuming that no 
other oxygen atom is at this position (see also Ref. [15]). Note that this type of 
motion is very different from the one found in the relaxation of a fragile glass-former 
(Lennard- Jones system [12]) since in that case the breaking up of the cage occurs 
in that the main peak observed in the (3— relaxation becomes broader and moves 
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FIGURE 4. Time dependence of the incoherent intermediate scattering function for the oxygen 
atoms for different temperatures. 

to larger distances and no secondary peak is observed at any time, thus showing 
that in these systems no jump processes are present. 

The relaxation dynamics of the system can also be studied well by means of the 
(incoherent) intermediate scattering function F"(q,t), which is given by a space 
Fourier-transform of Gf(r,t), i.e. 

i?(?,0 = ^« 1 E< ex p(^- (^(0-^(0)))) • (5) 

The time dependence of F°(q,t) is shown in Fig. 4 for all temperatures investi- 
gated. The value of q is 1.7A -1 , the location of the first peak in the static structure 
factor. We have found that the correlation functions for different wave-vectors as 
well as those for the silicon atoms are qualitatively similar to the one presented here. 
The same is true for the coherent intermediate scattering functions. From the figure 
we see that with decreasing temperature the a— relaxation time increases quickly 
and we have found that its temperature dependence is very similar to the one of 
the diffusion constants (see Fig. 2) in that it shows an Arrhenius dependence at 
low temperatures and at high and intermediate temperatures a dependence which 
can be fitted well with the power-law given in Eq. (3) [11]. (It was from this tem- 
perature dependence that the critical temperature T c = 3330 K was determined.) 
In particular we find that for temperatures above T c , the product between r and D 
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FIGURE 5. Incoherent intermediate scattering function for the oxygen atoms as a function of 
rescaled time t/r. 

was essentially constant, as it is predicted by MCT, whereas below T c the product 
increases (which is not in contradiction to the extended version of MCT). 

We see that also in F s (q,t) the cage effect is seen very well in that at low tem- 
peratures the relaxation function shows a plateau. Before the correlators reach 
this plateau, they show a pronounced dip, a feature not present in the dynamics 
of fragile glass-formers. This dip is the signature of the so-called boson-peak, a 
dynamical feature whose origin is still a matter of debate [9,16]. 

For the time-range of the a— relaxation, MCT predicts that the so-called time- 
temperature superposition principle (TTSP) holds, i.e. that a time correlation func- 
tion <p(t) can be written as 



where $ is a master function which will depend on 0. This prediction of the theory 
can easily be tested by plotting the time correlation function versus the rescaled 
time t/r, which is done in Fig. 5 for the intermediate scattering function shown 
in Fig. 4. From the figure we see that in the a— relaxation regime this type of 
scaling leads indeed to a nice collapse of the curves for the different temperatures, 
as predicted by the theory. In the time regime of the late j3— relaxation, however, 
strong deviations from the TTSP are observed, in contrast to what has been found 
in many fragile glass formers. A detailed analysis of the curves has shown [14] 



4>(t) = Ht/r(T)) 



(6) 



that the violation of the TTSP in this time regime is due to the presence of the 
boson-peak, since the intensity of this feature depends significantly on temperature 
and thus leads to a strong temperature dependence of the height of the plateau 
and hence to a violation of the TTSP in this time window. However, the shape of 
the curves was found to be essentially independent of temperature, in agreement 
with the prediction of MCT [14]. In particular we have found that in the late 
f3— relaxation regime this shape, i.e. the time dependence, is described very well by 
the so-called von Schweidler law, i.e. F s (q,t) = f%(q) — h(q)t b , where according to 
MCT the exponent b is related to the critical exponent 7 of Eq. (3) [1], a relationship 
which we have found to hold to within 10% for the present system. 

It is sometimes argued that the ideal version of MCT does not hold below T c , since 
real systems will always show ergodicity restoring hopping processes, and since the 
details of the extended theory have not been worked out, the whole theory becomes 
useless at low temperatures. In the following we will demonstrate that such an 
attitude is overly pessimistic and does not correspond to reality. Although the 
extended version of MCT does not make predictions on the a— relaxation regime, 
it makes useful and testable prediction in the /3— relaxation regime. In particular 
the theory predicts that in this regime the so-called factorization property should 
hold. This means that any time correlation function cj>{t) can be written as follows: 

(/>(t) = </> c + hg(t) . (7) 

Here the constant C is the height of the plateau (see Fig. 4), also sometimes called 
nonergodicity parameter, and h is a constant which is called critical amplitude. 
Both quantities depend on <fi. The important point is that the time dependence 
of (j) is given by the system universal function g(t), i.e. a function that does not 
depend on <fi. One possibility to test whether the factorization property holds is to 
calculate the ratio 

m) ~ W) - W) ' ( ) 

Here t' and t" 7^ t' are two arbitrary times in the {3— relaxation regime. A trivial 
calculation shows that the quantity R^{t) will be independent of if (pit) can be 
written in the form given by Eq. (7). In order to test this prediction of the theory 
we therefore plot in Fig. 6 the time dependence of R<j, for different correlators <p(t) 
at T = 2750 K. In particular we use for cf>{t) the incoherent intermediate scattering 
function for oxygen at q — 2.8A" 1 , the coherent intermediate scattering function 

for the Si-Si correlation at q = 1.7 A' 1 , the function b s (t) = Jq AA Anr 2 Gf(r,t)dr, 



and a Si _ (t) = ^r 2 G^M) - 1 



[0Si-o(r) - I] dr. Here Gf~°(r,t) and 
(?Si-o( r ) are the distinct part of the van Hove correlation function and the radial 
distribution function between Si-0 atoms, respectively [17]. (These two last func- 
tions have been proposed by Roux et al. and have no deep physical meaning, but 
are just functions which can be calculated with relatively high accuracy [15]. The 
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FIGURE 6. The ratio R^(t), see Eq. (8), for different correlation functions ^ as a function of 
time. T = 2750 K. See text for details. 

occurring integration boundaries are related to the location of the minima in the 
radial distribution function [14].) 

From the figure we see that in the time scale of the /3-relaxation regime the 
different curves collapse nicely onto a master curve, as it is predicted by MCT. For 
times outside this regime no collapse is observed, demonstrating that the existence 
of a master curve is a nontrivial feature. We emphasize that the curves are for a 
temperature which is considerably below T c and that, as demonstrated above by 
means of the diffusion constants and the van Hove correlation functions, at this 
temperature the hopping processes are very prominent. Thus we conclude that, 
although not all details of the extended version of MCT have been worked out, this 
theory is indeed able to make interesting and testable predictions also below T c . 
Finally we mention that we have found that the factorization property also holds 
well at temperatures as high as T = 4000 K, i.e. above the critical temperature 
T c = 3330 K, in agreement with the theory. 



We have presented the results of a large scale molecular dynamics computer sim- 
ulation of a realistic model, the BKS model [6], of supercooled silica, the prototype 
of a strong glass-former. Although we find that the time dependence of the mean 



SUMMARY AND DISCUSSION 



squared displacement of a tagged particle is roughly similar to the one in a fragile 
glass-former [12], important differences can be noted, such as, e.g., an overshoot 
which is related to the boson-peak as well as the temperature dependence of the 
long time behavior. The latter is related to the temperature dependence of the dif- 
fusion constant D and we find that D{T) shows at low temperatures an Arrhenius 
behavior with activation energies that agree very well with experimental values, 
thus giving evidence that this silica model is indeed quite realistic. With increas- 
ing temperature significant deviations from the Arrhenius-law are observed and we 
find that the temperature dependence of D can be fitted well with the power-law 
proposed by MCT with a critical temperature T c around 3330 K. Thus we find that 
for this strong glass former T c is more than twice the experimental glass transition 
temperature T g which is at 1450 K. (We note that in Ref. [11] we have given evi- 
dence that the BKS model has an experimental glass transition temperature which 
is close to the real one.) Thus these results support the idea [18] that also for 
strong glass-formers there is a temperature range in which the predictions of the 
ideal version of MCT hold, that, however, for the case of silica this range is at 
temperatures which are currently difficult to access by experiments. 

For temperatures below T c we find, by investigating the self part of the van 
Hove correlation function, that the motion of the particles is strongly influenced 
by jump like motions of the oxygen atoms, i.e. that instead of the flow-like motion 
of the particles proposed by the ideal version of MCT above T c , it might be more 
appropriate to think of a quasi-frozen potential energy landscape in which the 
dynamics is mainly given by the hopping of the particles [19]. That the presence of 
these strong hopping processes does not exclude the possibility that MCT makes 
useful predictions also at low temperatures, is demonstrated by our investigation of 
the factorization property [see Eq. (7)]. We show that this property holds very well 
even at temperatures at which the dynamics of the system is strongly influenced by 
hopping processes. Furthermore we have found [11] that the master curve obtained, 
see Fig. 6, is in the late part given by the von Schweidler law with an exponent b, 
which fulfills the relation proposed by MCT between b and the critical exponent 7 
of the relaxation times [1]. 

Finally we have also investigated the time and temperature dependence of the 
incoherent intermediate scattering function. We find that the time-temperature 
superposition principle holds well in the late a— relaxation regime that, however, 
quite large deviations are observed in the /3— relaxation regime. The reason for 
these differences is the boson-peak, which leads to a temperature dependence of 
the height of the plateau. The form of the curves, however, is indeed independent 
of temperature, as can be seen from the factorization property. 

Thus we can summarize by saying that many aspects of the dynamics of this 
strong glass-former can be understood by means of MCT, although the agreement 
between the ideal version of the theory and the results of the simulation is less 
satisfying as it has been found to be for fragile glass-formers. In particular the range 
over which the power-law for the diffusion constants can be observed is relatively 
small and the strong influence of the boson-peak onto the (3— relaxation leads to a 



breakdown of some, but not all, of the predictions of the theory in this time window. 
Nevertheless, the extended version of the theory is able to describe correctly some 
aspects of the dynamics also at low temperature, a result which should encourage 
to work out this version of the theory in more detail. 
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